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^} ■ Abstract: We show that the ABJM theory, which is an M = 6 superconformal 

U(iV) x U(iV) Chern-Simons gauge theory, can be studied for arbitrary N at arbitrary 
coupling constant by applying a simple Monte Carlo method to the matrix model that can 
be derived from the theory by using the localization technique. This opens up the possi- 
bility of probing the quantum aspects of M-theory and testing the A0IS4/CFT3 duality at 
the quantum level. Here we calculate the free energy, and confirm the iV 3 / 2 scaling in the 
M-theory limit predicted from the gravity side. We also find that our results nicely inter- 
polate the analytical formulae proposed previously in the M-theory and type IIA regimes. 
Furthermore, we show that some results obtained by the Fermi gas approach can be clearly 
understood from the constant map contribution obtained by the genus expansion. The 
method can be easily generalized to the calculations of BPS operators and to other theo- 
ries that reduce to matrix models.* 
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1. Introduction 

Three-dimensional gauge theories with a Chern-Simons term have been studied extensively 
for their rich and beautiful mathematical structure [1] and for applications to quantum 
Hall systems in condensed matter physics. Their relevance to superstring/M-theory was 
realized in 2004 [2] by the appreciation that superconformal Chern-Simons theories coupled 
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to matter fields capture the dynamics of M-theory on a multiple M2-brane background. The 
latter theory is expected to be obtained as an infrared fixed point of type IIA superstring 
theory on a D2-brane background, whose low-energy description is given by the maximally 
supersymmetric (2 + l)-dimensional Yang-Mills theory. In 2008, an explicit form of the 
Chern-Simons theory that describes the infrared fixed point was proposed by Aharony, 
Bergman, Jafferis and Maldacena (ABJM) [3], following important earlier works [4, 5]. It 
is a U(N) x U(N) theory with Chern-Simons levels k and —k coupled to bifundamental 
matters. The on-shell supersymmetric Lagrangian of the theory is given by 

£-U(N)xU(N) 



fermionic complex bifundamental fields, respectively. This theory has N = 8 supersymme- 
try for k = 1,2 and N = 6 supersymmetry for k > 3. It has been conjectured to be dual 
to M-theory on AdS± x S 7 /Z k for k < iV 1 / 5 , and to type IIA superstring on AdS± x CP 3 
in the planar large- N limit with the 't Hooft coupling constant A = N/k kept fixed. 

From the viewpoint of quantum gravity, the ABJM theory is important since it may 
provide us with a nonperturbative definition of type IIA superstring theory or M-theory 
on AdS± backgrounds since the theory is well-defined for finite N. One may draw a pre- 
cise analogy with the way maximally supersymmetric Yang-Mills theories may provide us 
with nonperturbative formulations of type IIA/IIB superstring theories on D-brane back- 
grounds through the gauge/gravity duality [6, 7, 8, 9]. In particular, the M-theory limit 
is important given that M-theory is not defined even perturbatively, although there is a 
well-known conjecture on its nonperturbative formulation in the infinite momentum frame 
in terms of matrix quantum mechanics [10]. The planar limit, which corresponds to type 
IIA superstring theory, has interest on its own since it may allow us to perform more de- 
tailed tests of the gauge/gravity duality than in the case of AdS^/CFT^. In particular, we 
may hope to calculate the 1/N corrections to the planar limit, which enables us to test the 
gauge/gravity duality at the quantum string level, little of which is known so far. 

In all these prospectives, one needs to study the ABJM theory in the strong coupling 
regime. As in the case of QCD, it would be nice if one could study the ABJM theory on a 
lattice by Monte Carlo methods. This seems quite difficult, though, for the following three 



k Tr p ( A^dyAp ^ApA v Ap + Apd v Ap + ^ApA u Ap^j 




(1.1) 



where Ap and Ap are U(N) gauge fields, and $ Q and ^ a (a = 1,2,3,4) are bosonic and 
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reasons. Firstly, the construction of the Chern-Simons term on the lattice is not straight- 
forward, although there is a proposal [11, 12] based on its connection to the parity anomaly. 
Secondly, the Chern-Simons term is purely imaginary in the Euclidean formulation, which 
causes a technical problem known as the sign problem when one tries to apply the idea of 
importance sampling. Thirdly, the lattice discretization necessarily breaks super symmetry, 
and one needs to restore it in the continuum limit by fine-tuning the coupling constants of 
the supersymmetry breaking relevant operators. (See, for instance, ref. [13].) This might, 
however, be overcome by the use of a non-lattice regularization of the ABJM theory [14] 
based on the large- ./V reduction on S 13 [15, 16], which is shown to be useful in studying the 
planar limit of the 4d M = 4 super Yang- Mills theory [17, 18, 19]. 

What we do here instead is to apply Monte Carlo methods not to the original theory 
(1.1) but to a matrix model obtained after a huge reduction of the degrees of freedom 
due to supersymmetry. In fact, it has been known for a while in certain super symmetric 
theories that one can reduce the path integral to a finite dimensional matrix model by 
using the so-called localization technique. Such a technique was applied [20] to 4d N = 4 
super Yang-Mills theory, and some conjecture on the half-BPS Wilson loops 1 [21, 22] has 
been confirmed. In ref. [23], the same technique has been applied to the ABJM theory on 
three-sphere S 3 , and its partition function was shown to reduce to a matrix integral 

1 r d N a d N u 
Z(N,k) - [ a ii a v 



(N\) 2 J (2it) n (2it) n 

n i<J (2sinh^) 2 (2sinh. 



2 cosh • 



exp 



§f>. 2 -^ 



i=l 



(1.2) 



which is commonly referred to as the ABJM matrix model. 2 By using this matrix model, 
the free energy of the ABJM theory has been studied intensively [28, 29, 30, 31, 32, 33, 
34, 35]. In ref. [29] the planar limit and the 1/N corrections 3 around it have been studied 
employing a technique from topological string theory, and the on-shell action of the type 
IIA supergravity on AdS^ x CP 3 has been reproduced. In ref. [30] the free energy in the M- 
theory limit has been obtained using some ansatz for the eigenvalue distribution. In ref. [33] 
the genus expansion at strong 't Hooft coupling has been considered and a resummed form 
was obtained in terms of the Airy function by using the holomorphic anomaly equation 
[36]. The obtained simple form was claimed to be valid to all orders in the genus expansion 
up to the worldsheet instanton effect. In ref. [35], the free energy in the M-theory regime 
at small k has been calculated by the Fermi gas approach, and the result turns out to be 
given by the Airy function obtained in ref. [33] with some extra terms. These results, if 



1 This formula is also reproduced by a numerical simulation in the large- N limit [17, 18]. 

2 The localization of the ABJM theory and related theories on various manifolds such as S s [24], S 1 x S 2 
[25] and squashed S 3 [26] has also been considered. 

3 There is also a study of non-planar corrections to anomalous dimensions by the integrability approach. 
See e.g., [27]. 
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correct, would enable us to shed light on the dynamical aspects of M-theory and to test 
the AdS/CFT duality including the string loop effect by studying the gravity side further. 

In this paper we show that the ABJM matrix model can be rewritten in a form suitable 
for Monte Carlo simulations, which enables simple calculation of the partition function and 
BPS operators for arbitrary values of the rank N and the level k from first principles. In 
particular, we calculate the partition function explicitly for various N and k, which is sup- 
posed to contain the nonperturbative effects corresponding to the worldsheet instantons in 
string theory neglected in refs. [29, 33]. We find the well-known constant map contribution 
[36, 37, 38] is also correctly reproduced. The constant map contributions are of the order 
of X° /k 2d ~ 2 , where g is the genus. Note that these terms depend on the string coupling 4 
g s = 2iri/k but not on the 't Hooft coupling constant. In the planar limit, they correspond 
to a constant shift of the free energy, which becomes negligibly small in the large 't Hooft 
coupling limit, and hence they do not spoil the agreement with the type IIA supergravity. 
They are negligible in the eleven-dimensional supergravity limit as well. However, when 
one considers quantum strings or M-theory, these terms should be taken into account since 
they have an explicit g s dependence. In fact, these terms become dominant at strong 't 
Hooft coupling for g > 2. 

This paper is organized as follows. In section 2 we review the previous results for the 
free energy of the ABJM matrix model in various limits, which are obtained by analytical 
methods. In section 3 we describe our numerical method. In section 4 we present our 
results, and discuss the discrepancies from the analytical results. In section 5 we show 
that these discrepancies can be interpreted as the constant map contributions. Section 
6 is devoted to a summary and discussions. In Appendix A we explain the basics and 
some details of Monte Carlo simulation. In Appendix B we derive an alternative form of 
the ABJM matrix model, which is suitable for Monte Carlo simulation. In appendix C 
we show the equivalence between the constant map contribution and the Fermi gas result 
A{k) — ilog2, which is derived by the large-/c and small-A; expansions, respectively. 

Note added. After the first version appeared on the arXiv, we were informed by Marcos 
Marino that the discrepancies observed at genus 0, 1 and 2 between our numerical results 
and the formula proposed in ref. [33] can be naturally attributed to the constant map 
contributions. We greatly appreciate this important comment, which enabled us to deepen 
our argument in section 5.2. 

2. Previous analytical results for the free energy 

In this section we summarize some known analytical results for free energy of the ABJM 
theory, which is defined in terms of the partition function (1.2) as 

F{N,k) = \ogZ{N,k). (2.1) 

4 This is the string coupling in the context of topological string theory. The string coupling in the IIA 
superstring limit of the ABJM theory is g a IIA ~ V\/k 2 ~ VXg 2 s . 
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2.1 Perturbative results for all N 



The free energy can be calculated by using a usual perturbative technique, and the result 
at the one-loop level is given as (See, for example, ref. [32].) 



2N 

Fweak = -iV 2 log— -iVlog2vr + 21ogG 2 (Ar + l) 

TTA 

JV>1 



(2.2) 



V ^ iV 2 (log 2vrA - I - 2 log 2 ) - I log N + 2 C '(-1) + £ ^T^)^ ^ 

g=2 

where G 2 (x) is the Barnes G-function G 2 (x) = n«=i s - The 1/iV-expansion is shown in the 
second line with Riemann's zeta function (,(x) and the Bernoulli numbers Big. The 0(N 2 ) 
terms in (2.3) agree with the result (2.5) obtained in the planar limit. Note, however, that 
the expression (2.2) includes contributions to all orders in the 1/iV-expansion. 

2.2 N = 2 with arbitrary k 

An exact expression for N = 2 is obtained by Okuyama [34] as 5 



F(2,k) = { 
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4 log 2 



4 log 2 for odd k 



(2.4) 



for even k . 



This result has been obtained by direct integration of (B.4). Since the expressions 
for the odd and even k cases are different, the analyticity in k (when one regards k or 
equivalently the 't Hooft coupling constant A as a continuous variable) is not obvious a 
priori. However, as we will see in section 4, our numerical results suggest that the free 
energy is a smooth function of k. The analyticity is important in the context of the 
AdS/CFT correspondence, in which one assumes the analyticity on the gravity side. Also 
the analysis in the planar limit usually assumes the analyticity implicitly. 

2.3 Planar limit (N ->• 00 with A fixed) 

The free energy in the planar limit (N — > 00 with A fixed) has been calculated by Drukker, 
Marino and Putrov (DMP) [29]. These results have been obtained by a standard matrix 
model technique after the analytic continuation [28] to the lens space L(2, 1) = S* 3 /Z2 
matrix model [39, 40], which is obtained from the pure Chern-Simons theory on L(2, 1). 
The validity of the analytic continuation is proved diagramatically in refs. [41, 42]. 
At weak coupling (A <C 1) the authors obtain 



r;,,k. ; .h,,,.., - N 2 (log2^A- I - 2 log 2 



(2.5) 



3 Note that the normalization of the partition function adopted in ref. [34] differs from OUrS aS ^Okuyama — 
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up to 0(A). 

At strong coupling (A> 1) the authors obtain 

7T\/2A 3/2 1 

^dmp = 1 -^N 2 where A = A - — (2.6) 

to all orders of the 1/A expansion. The leading behavior Fdmp — — v / 2vriV 2 /(3v / A) agrees 
with the dual type IIA supergravity prediction [29, 43] including the overall coefficient. It 
has been claimed that the free energy (2.6) at strong coupling receives the correction of 
the form 

A 7m \W2\J 

where ff\x) is a polynomial in x of degree 21 — 3 (for I > 2). This exponentially small 
correction has been interpreted in ref. [29] as the effect of the worldsheet instanton in the 
dual type IIA superstring, which corresponds to a string worldsheet wrapping a CP 1 cycle 
in CP 3 [45]. 

In section 4 we will show that another contribution of the order of 0(N 2 /X 2 ) due to the 
constant map needs to be added in comparing with precise numerical analysis. Although 
this term does not affect the agreement with supergravity, it must be taken into account 
when one compares the finite A corrections with the string a' corrections. 

2.4 M-theory limit (N ->• oo with k fixed) 

In ref. [30], the free energy in the M-theory limit 6 (N — > oo with k fixed) has been calculated 
and confirmed the prediction 



FsuGRA = -^iV 3 / 2 (2.7) 

from the dual eleven-dimensional supergravity, which shows the well-known N 3 ^ 2 scaling 
for the degrees of freedom in the theory of M2-branes [44]. Note also that (2.7) agrees with 
what one obtains formally from the leading large-A behavior of the planar result (2.6) by 
replacing A with N/k. 

The result (2.7) was obtained by imposing an ansatz for the eigenvalue distribution 

Hi = N a Zi + iwi , Vi = N a Zi - iwi (zi,Wi £ R) , 

which is necessary for the cancellation of long-range forces, and is also suggested by numer- 
ical studies of the saddle point equation. The parameter a is chosen to be 1/2 by requiring 
that all the short-range forces contribute to the free energy at the same order of N in order 
to have nontrivial solutions. 



6 Strictly speaking, since the ABJM theory has been conjectured to be dual to the M-theory for k <S N 1 ^ 5 , 
the limit N — > oo with k fixed is merely a sufficient condition. In the following, however, we simply call it 
"the M-theory limit" . 
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2.5 1/N expansion around the planar limit 

Fuji, Hirano and Moriyama (FHM) [33] studied the free energy to all orders in the genus 
expansion neglecting the instanton contribution, which is of the order of 0(e~ 2?rv/ ^). Their 
proposal for a resummed form is given by 

2/3" 



Ffum(N,X) = log 



J_ / 47T 2 iV \ 1/3 



Ai 




(2.8) 



where Ai(x) is the Airy function, and the "renormalized 't Hooft coupling" A rcn is given by 

A rcn = A-l-^. (2.9) 

The appearance of the Airy function [33] is also encountered in the context of M-theory 
flux compactification [46]. Note that the expression (2.8) reproduces (2.6) in the large- 
N limit as one can easily see by using the asymptotic formula logAi(x) ~ — 2x ^ /2 for 
i>1, In section 4 we will show that (2.8) has another contribution, which is necessary 
for comparison with our numerical results. 

The free energy at higher genus has been studied earlier [29, 31] by using a topological 
string technique after analytic continuation to the lens space matrix model. The analysis in 
ref . [33] has been performed by using the holomorphic anomaly equation [36] , whose solution 
is the same as the one for the loop equation [47, 48] with some appropriate boundary 
conditions. In order to solve the holomorphic anomaly equation, one needs to provide 
some inputs such as the free energy at genus zero and one, which are taken to be 

fffi, = *f=V and F& l = ^V*-\v*(&). 

In this way the authors have found a general solution, which gives the free energy at all 
genus up to the worldsheet instanton effect. The integration constants were determined by 
assuming the absence of non-perturbative corrections of the type ~ O Strictly 
speaking, what one obtains in this way is the "weight zero" contribution to the free energy 
in the language of topological string theory. It is claimed that one can turn this result into 
the one including contributions from all weights by making a replacement A — > A rcn , which 
is given in (2.9). 

This "renormalized 't Hooft coupling" is different from the expectation from the gravity 
side [49] : A ren ,grav = A - 1/24 + A 2 /(24iV 2 ). While it is possible that this disagreement 
may imply that the AdS/CFT does not hold at finite- iV/quantum string level, we should 
definitely gain more understanding on both gauge theory and gravity sides. The additional 
contribution to the FHM result from the constant map should be important also from this 
point of view. 

2.6 N » 1, small k 

In ref. [35], the free energy with fixed small k has been calculated by using the Fermi gas 
approach neglecting the quantum mechanical instanton effect (worldsheet instanton) and 
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the terms which are suppressed exponentially at large iV (membrane instanton). In this 
approach, the partition function of the ABJM theory is regarded as an ideal Fermi gas 
system described by (B.2) with the Planck constant identified as h = 2irk. The result is 
given by 

(47T 2 fc) 1/3 [/7rfc 2 \ 2/3 (N 1 J_ 

^2 1 \V2J \k 24 3k 2 



^Fe 



log 



+ A(fc)--log2. (2.10) 



The leading large- N behavior reproduces eq. (2.7) exactly. The function A(k) in (2.10) is 
given for k <C l/(27r) as 7 



A(k) = 



2C(3) 



k_ 
12 



7T 2 fe 3 



0(k 5 



(2.11) 



7r 2 k 12 4320 

Since the first term in (2.10) can be obtained formally from the FHM result -Ffhm in 
(2.8) by replacing A with N/k, one can rewrite it as 



^Fcrmi = i*FHM + A(k) - - log 2 . 



(2.12) 



where A{k) may be viewed as "quantum corrections" with the "Planck constant" h = 2nk. 
Note that the first term in (2.12) is valid for all k although (2.11) is obtained at small k. 
The authors note that the second and third terms in (2.11) are given by 



A(k) = 



203) 

TT 2 k 

203) 
ir 2 k 



E 



2n 



' n(2n 



-IT 



2n-2^,2n-l 







l)(2n)! 

sin(e/2) 



(2.13) 



for n = 1,2. This is the power series with odd powers of k unlike the usual genus expansion 
around the planar limit. The authors suggest that A{k) may encode the effect from D0- 
branes of the order of O (e~ k ) ~ O (e -1 / fls ). 

Since this analysis for A(k) assumes k <C l/(27r), it is not clear a priori whether the 
result holds at physical values of k corresponding to integers. As we will see later, (2.11) and 
(2.13) are in reasonable agreement with our numerical result for small k such as k = 1, 2, 3, 
but not for larger k (including the planar limit). 



3. Numerical methods for the ABJM matrix model at arbitrary N and k 

In this section we discuss how we can study the ABJM matrix model at arbitrary and 
k by applying a standard Monte Carlo method. For the readers who are not familiar with 
Monte Carlo methods in general, we review the basic ideas in Appendix A. For an earlier 
work on Monte Carlo simulation of a one-matrix model, see ref. [50]. 

7 Although the Chern-Simons level k must be integer in a physical setup, the integral (1.2) is itself well- 
defined also for non-integer k and we can actually obtain numerical results, which turn out to be a smooth 
function of k. 
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The ABJM matrix model in the form (1.2) is not suitable for Monte Carlo simulation 
since the integrand is not real positive. 8 However, as we review in Appendix B in detail, 
one can rewrite the ABJM matrix model as follows. 

Z(N,k) = C Ntk g{N,k), C NM = 



g(N,k) = d N x^ } 2 \ 7 . (3.1) 

K ' J IL 2 00811(^/2) V ; 

In the k = 1 case, one may view (3.1) as a mirror description of the ABJM theory in 
terms of the 3d U(N) N = 4 SYM with adjoint and fundamental hypermultiplets, which 
is isomorphic to 3d U(N) N = 8 SYM in the low-energy limit [51]. The important point 
here is that, in this form (3.1), the integrand is real positive, and we can perform Monte 
Carlo simulation in a straightforward manner as described in Appendix A. 

We should also note that, while the level k should be an integer in the original 3d 
gauge theory, nothing prevents us from considering non-integer k in the integral (1.2). In 
what follows, we therefore extend the value of k to any real number. 

In order to calculate the free energy (2.1), which is the log of the partition function, 
we need to rewrite it in terms of expectation values of some quantities, which are directly 
calculable by Monte Carlo methods. 9 The basic idea in our case is to calculate the ratios 
of the partition functions for different k or N as expectation values. Since we know the 
results for k = or TV = 1, we can obtain results for arbitrary k and N by calculating 
an appropriate product of the ratios. Depending on whether we change k or N, we have 
the following two methods, which give the same result within statistical errors as we have 
checked for various k and N. The second method is particularly useful in studying the 
M-theory limit, which corresponds to the large N limit with fixed k. 

3.1 Calculating the ratio of partition functions with different k 

Let us consider a trivial identity 

9(NM) = fd N xe- S ^) / —S(N,k2;x)+S(N,ki;x)\ ,3 2) 

g(N, fei) jd N xe-S{NM;x) \ e / NM ' ^ > 

where we have defined 

-S(N,k;x) = Ui<j ta:ah2 i. ?h w 2 -) , s 

n i 2cosh(x J /2) 1 • ) 

and ( • • • ) jv.fc stands for the expectation value with respect to the action S(N, k; x) 

<°>".* = J d N xe -S(N,k;x) ■ ( 3 - 4 ) 



8 0ne might think of simulating a system without the phase factor exp((ifc/47r)(/4 — uf)), and including 
its effect afterwards by reweighting. While it is possible to obtain results for the k = 1 case along this line, 
the calculation becomes more and more difficult for larger k due to the sign problem. 

9 For applications of such an idea on different supersymmetric systems, see refs. [52] and [53]. 
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The quantity (3.2) can be calculated easily by the standard Monte Carlo method as far as 
k\ and k 2 are sufficiently close. 10 Therefore, we can calculate the free energy F as 

F = log Z = log CV,n + log 9 (JV, k) 

= l„gO N , l + gl„g-|^L + ,og 9W 0) 

I 

= log C N , k + V log / e -SW*i?)+S(Af-*i-ii*)\ + TV log vr , (3.5) 

where = ko < ki < ■ ■ ■ < ki = k and we have used g(N, 0) = J r~[ 2cosh(a: /2) = ^ in tne 
last line. We have to make the adjacent values of k close enough for the reason mentioned 
above. 

3.2 Calculating the ratio of partition functions with different N 

Let us decompose N into N = N\ + N2 and consider the ratio 

g(N,k) _ j d jV xe -S(iV,fc) 

g(N 1 ,k)g(N 2 ,k) ~ J d N x e -S(JVi,fc;xi,- ,x Nl )-S(N 2 ,k;x Nl+1 ,- ,x N ) 

— /gS(JVi,fc;xi,- ,3j JVl )+5'(A r 2,fc;a;jv 1 +i,'-- ,x N )-S(N,k)\ (3 6) 

where the symbol {•••)n 1 ,n 2 denotes the expectation value with respect to the "action" 
S(Ni, k; xi, ■ ■ ■ ,XN 1 ) + S(N2,k;xN 1 +i,-" , £jv)- Note that 

JVi JV / N 
e S(N u k;xu-,x Nl )+S(N 2 ,k;x Nl + 1 ,-,x N )-S(N,k) = "Q "Q tanh 2 / ^ ^7 \ ^ 

due to the factorization of the potential terms. In order to calculate the right-hand side of 
(3.6) with good accuracy, it is necessary to take N2 small enough to make sure that (3.7) 
does not fluctuate violently during the simulation. In actual calculation we use N2 = 1. 
Then by calculating (3.6) for Ni = 1, 2, 3, ■ ■ ■ , and by using the N = 1 result 

tw-JiJm-'- (3 - 8) 

we can calculate the free energy for N = 2, 3,4, • • • successively with a fixed value of k. 
4. Results for the free energy 

In this section we present our numerical result for the free energy of the ABJM theory. 
In order to test our code, we first study the N = 2 case and compare our result against 
the exact result (2.4) obtained by Okuyama [34]. As can be seen from fig. 1, our result 
reproduces the exact result very accurately. We have also obtained results for non-integer 
values of k, which are not obtained in ref. [34]. They are found to connect the results for 
integer k smoothly. 

10 As k 2 moves away from fei, the quantity e -S( Jv .' ! 2;a:)+s(JV,fci;a:) fj uc tuate violently during the simulation 
of the system S(N, ki;x), which leads to larger statistical errors. 
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Figure 1: The free energy of the ABJM theory for N = 2 is plotted against the Chern-Simons level 
k. The circles and triangles represent our Monte Carlo result and the exact result (2.4), respectively. 

4.1 Planar limit 

Next we consider the planar limit (N — > oo with A = N/k fixed), which is conjectured 
to be dual to the classical type IIA superstring on AdS^ x CP 3 . In fig. 2 we plot the 
normalized free energy F/N 2 against 1/N 2 for various values of A. Our results can be 
fitted well by F(N, \)/N 2 = / (A) + fi(X)/N 2 - ± log N as theoretically expected 11 . In the 
left panel of fig. 3, we plot /o(A) = lim j /v_ s>00 F(N, \)/N 2 against 1/VX The results seem 
to interpolate the DMP result (2.6) at strong coupling and the perturbative result (2.5) at 
weak coupling. However, by looking more carefully into the asymptotic behavior for large 
A, we find certain discrepancies. This can be seen from the right panel of fig. 3, in which 
we plot the difference limjv_ s . 00 (i ? — -Fdmp)/-^ 2 ! which is found to behave as 

F — F DM p a»i a 
ife N 2 " A^ + 6 °' (41) 
a = -0.015 ± 0.001 , b = -0.0006 ± 0.0002 (4.2) 

instead of the behavior 0(e^ 2nv ^) expected from the worldsheet instanton effect. We 
consider that bo is consistent with zero since the fitting error may well be slightly underes- 
timated. Since the discrepancy (4.1) vanishes at A = oo (assuming that bo in (4.1) is zero), 
it does not affect the agreement with the dual type IIA supergravity. 

In section 5 we explain that this discrepancy can be understood as the constant map 
at genus 0. Similar discrepancies exist also in 1/N corrections around the planar limit as 
we will see. 

4.2 M-theory limit 

Next we consider the large- A" limit with fixed k, which is conjectured to correspond to the 
eleven dimensional supergravity on AdS^ x S"' ' /^k- Figure 4 shows that the free energy 



ii 



as 



The functions /o(A) and /i(A) defined here are related to Fo(\) and -Fi(A), which are defined in (5.1), 
/o(A) = -F (A)/4tt 2 A 2 and /i(A) = Fi(A). 
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Figure 2: The normalized free energy F/N 2 is plotted against l/N 2 for various values of A 
(Left). In the right panel, we zoom up the plot for A — 1. The data can be nicely fitted to 
F(N, X)/N 2 = ./o(A) + fi(X)/N 2 — i log TV, which enables us to make a reliable extrapolation to 
the planar N — > oo limit. 




Figure 3: (Left) The free energy in the planar limit Jo (A) = lhriAr^oo F(N, X)/N 2 extracted from 
fig. 2 is plotted against 1 / \f\. Our results seem to interpolate the DMP result at strong coupling 
and the perturbative result at weak coupling. (Right) The difference between our result and the 
DMP result, i.e., lim7v->oo(-F 1 — -Fdmp)/TV 2 , is plotted against 1/A 2 . The data points can be fitted 
to a straight line, which implies (4.1) and (4.2). 

b grows in mag nitude as iV 3/2 with N, and F/N 3 / 2 behaves as F(N,k)/N 3 / 2 = h (k) + 
h\{k)/N, which enables us to obtain the M-theory limit ho(k) = limAr_ >00 F(N, k)/N 3 / 2 
reliably. 

In fig. 5 we plot ho(k) against \^k, which confirms the prediction (2.7) from eleven- 
dimensional supergravity for k = 1, 2, • • • ,10 very precisely. 

4.3 Finite N effects 

One of the important results on finite iV effects in the free energy is that the l/N corrections 
around the planar limit are resummed in a closed form (2.8) neglecting the worldsheet 
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Figure 4: (Left) The free energy is plotted against N 3 / 2 for k = 1,2,4,6,8. The data points can 
be fitted to straight lines, which implies F ~ N 3 / 2 as N increases. (Right) The normalized free 
energy F/N 3 / 2 is plotted against l/N. The data can be nicely fitted to straight lines, which enables 
us to make extrapolations to the M-theory limit reliably. 
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Figure 5: The M-theory limit of the free energy limjv_ i . 00 F/N 3 / 2 extracted from fig. 4 (Right) is 
plotted against \fk. Our data are in good agreement with the result (2.7) predicted from eleven- 
dimensional supergravity, which is represented by the solid line. 

instanton effect. In fig. 6 we plot our results for N = 4 and N = 8 and compare them 
with the FHM result (2.8) and the DMP result (2.6). We find that both FHM and DMP 
are close to our data at strong coupling, but the difference between them is too small to 
see whether FHM is doing any better than DMP. This is simply because the term (2.7), 
which commonly exists in both results, is dominating over the difference. We therefore 
plot F — -Fsugra against N for k = 1 (Left) and k = 8 (Right) in fig. 7. The leading large- 
TV behavior of the plotted quantity is ^(53 + ^)A; 3/2 \/iV for FHM and -^/^k 3/2 VN 
for DMP, where the difference comes from the A 2 /(3iV 2 ) = l/(3k 2 ) term in (2.9). The 
difference becomes negligible for k = 8, but it is significant for k = 1, in which case our 
data are indeed closer to FHM than to DMP. 

We also find some discrepancy between our result and FHM, which are almost inde- 
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Figure 6: The free energy of the ABJM theory for TV = 4 (Left) and N = 8 (Right) is plotted 
against 1/vA. The solid line and the dashed line represent the FHM result (2.8) and the DMP 
result (2.6), respectively. The dotted line represent the perturbative results (2.2). 




2 4 6 8 10 12 14 16 2 4 6 8 10 12 14 16 

N N 

Figure 7: The free energy of the ABJM theory after subtraction of the dominant term (2.7) is 
plotted against TV for k — 1 (Left) and k — 8 (Right). The solid line and the dashed line represent 
the FHM result (2.8) and the DMP result (2.6), respectively. 



pendent of N. To see it more directly, we plot in fig. 8 the difference between our result 
and the FHM result for various k. It turns out that the discrepancies are indeed almost 
independent of N. This strongly suggests that the FHM result correctly incorporates the 
finite N effects except for a term which depends only on k. Note that this discrepancy 
cannot be explained by the worldsheet instanton effect 0(e~ 2n ^), which is neglected in 
FHM. While this discrepancy does not affect the M-theory limit corresponding to the strict 
N — > oo limit for fixed k, it is non-negligible when one considers 1/N corrections. As we 
will see in section 5, this discrepancy coincides with A(k) — ^log2 in eq. (2.12) by Fermi 
gas approach [35] for small k and with the constant map contribution for all k. 
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Figure 8: The difference F — -Ffhm is plotted against N for various values of k. It reveals non- 
negligible discrepancies for each k, which are almost independent of N. 

5. Interpretation of the discrepancies 

In this section we provide an interpretation of the discrepancies between our data and the 
known analytical results, which we observe in the previous section. 

5.1 Genus expansion 

Let us consider the planar limit, in which g s N = 2iriN ' jk = 2iri\ is kept fixed. In that 
limit, the free energy can be expanded with respect to the genus as 

oo 

F(g s ,X) = J2F 9 (X)9?- 2 

9=0 



V F (A) + Fi(A)-£^F a (A) + .... (5.1) 



(2ttX) 2 UV ' 1V J iV 2 
Below we consider the free energy order by order in this expansion. 

Planar contribution 

The planar contribution — fc 2 i ? o(A)/(4-7r 2 ) can be studied by the saddle point method, and 
the Fq(X) can be determined by solving [28, 29, 32, 54] 



cWA) = 7 G$ 



/ 1 1 1 

K n 2,3 2 2 2 



4 3,3 I -i 



K 2 \ i7T 2 K (I 1 1 3 K 2 . 

- +^- 3 F 2 (-,-,-;!,-;--) , (5.2) 



K „ (\ 1 1 3 K 2 



AW = ^ U'2'2 :1 -2 : -16.J ' (5 ' 3) 
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where G33 is the Meijer G-function 12 and 3F2 is the hypergeometric function. Note that 
these equations are exact for arbitrary A, and hence they fully incorporate the worldsheet 
instanton effect. One can obtain -Fo(A) by carrying out the integration over A as 

r\ /•k(a) a\i 

F (A) = F o (0) + J dX' dyF (X') = F (0) + J dn' — dyF (X') . (5.4) 

At weak coupling A <C 1, in particular, one obtains 

F (A) = F (0)+F (A) + O(A 9 ), (5.5) 

~ , 2 o/ vrA\ 4tt 4 A 4 61vr 6 A 6 12289vr 8 A 8 , , 

F„(A)=2^(3-21og T ) + - — + ^ m -. (5.6) 

By comparing this with the perturbative calculation (2.5), one finds Fq(0) = 0. 
At strong coupling A> 1, one obtains 



(5.7) 




(5.8) 

where A = A — 1/24. Here Co is an "integration constant", which has been set zero in the 
previous works, for instance, in ref. [29], which leads to eq. (2.6). 

In fig. 9 we plot c(A) = Fq(X) — Fq(X), where Fq(X) is evaluated numerically by per- 
forming the integral (5.4) explicitly. As A increases, we find that c(A) approaches a nonzero 
constant 

c - 0.60103, (5.9) 

which coincides with £(3)/2 ~ 0.601028 obtained as the constant map contribution at 
genus zero [36, 37, 38]. The value of ao in (4.1) predicted from the above calculation is 
a = -co/Ait 2 ~ -0.015224, which agrees well with the discrepancy (4.2) observed in the 
planar limit. 

Higher genus contributions 

Next we discuss the discrepancy at higher genus, which can be also interpreted as the 
constant map contributions as in the planar part. Let us note that the analytical results 



2 The Meijer G-function is defined by 

nr=i r fe- s )n; =1 r(i-a J + s 



( ai • 


• a p 








■) 




■ b q 





_L / : 

2m J L W j=m+1 r(i - h + 8 ) UU+i T ^ - s ) 



where the path L is chosen in an appropriate way depending on the parameters. See, for instance, ref. [55] 
for the details. 
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Figure 9: The solid line represents c(A) = Fq(X) — Fq(X), where Fq(X) is evaluated numerically 
by performing the integral (5.4) explicitly. The "integration constant" Co in (5.7) is obtained as 
(5.9) from the asymptotic value of c(A) at large A. The dotted line represents Fq(X) — Fq(A), where 
F (X) is the result at weak coupling given by (5.6). The matching of the weak coupling result and 
the strong coupling result around A ~ 0.15 also requires a similar value of Co in (5.7). 

up to the constant map have been obtained for genus one and two in terms of modular 
forms as [29, 31, 48, 56] 

^duiarW—logrKr), (5.10) 
F — r (A) = + 3tf^f - 2E 4 E 2 ) 



16^ 2 + 15g|gj + 2W%dj + 2$f 
12960^i9| 



(5.11) 



where rj(r) is the Dedekind eta function, E n (r) is the Eisenstein series of weight n, i? n (r) 
is the theta function, and t(A) is defined as 

= i K{f) = 1 + I^o(A) , (5.12) 

where K(x) is the complete elliptic integral of the first kind. 

In fig. 10 we plot F^ T - [f^ + j log k) (Left) and F^ dular - F^ (Right) 
against A, where we have defined the weak coupling results 

F £LkW = -^(logA + logfc) + 2C , (-l) , (5.13) 

at genus one and two, respectively, as can be read off from (2.3). We find in both cases 
that the result approaches a constant as A — > 0, which gives 

ci ~ 0.25558 , c 2 ~ 0.0027777 , (5.15) 
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Figure 10: (Left) The solid line represents f^dular — f-F^ak + I 1°§ ^) ' which is the difference 
between the modular expression (5.10) and the weak coupling result for genus one. It approaches 
a constant smoothly for A — > 0, which gives ci in (5.15). The dotted line represents -PpHM — 
(-^weik + \ 1°S ' w hich diverges as A — > since the FHM result neglects the worldsheet instanton 

(2) (2) 

effect. (Right) The solid line represents F 1 ^ > dular — -F^eak' which is the difference between the 
modular expression (5.11) and the weak coupling result for genus two. It approaches a constant 



smoothly for A — > 0, which gives C2 in (5.15). The dotted line represents i*pjjM — ^weak» which 

(2) 

diverges as A — > since the FHM result neglects the worldsheet instanton effect. Here Fp^ M is 



r(2) 



given by if H > M 



144^^ 



A-i/2 



1 x- 1 



A- 3 / 2 . 



respectively. This suggests that in the weak coupling region there are additional terms 
given by 



-C2 



AF (1) = — logfc-ci, AF (2) 
6 

which precisely coincide with the constant map contributions [36, 37, 38] 
for genus 1 : 
for genus g > 2 : 



(5.16) 



i log fc + 2^-1) + £ log-, 
6 6 2 

4 9 B2gB 2 g-2 



4g(2g-2)(2g-2)\ 



(5.17) 



Since the FHM result (2.8) reproduces the previous results in the genus expansion, the 
FHM result must also have the additional contributions 



k 2 1 47T 2 



(5.18) 



where the worldsheet instanton effect is neglected. 

As we did in the case of planar contribution, we can test whether the discrepancy 
in the genus one contribution between our data and the previous analytical results can 
be explained by the additional terms identified above. In fig. 11 we extract the genus 
one contribution from our data in the following way. First we subtract from our data for 
the free energy, the planar contribution <?~ 2 .Fo(A), where Fq(X) is obtained by (5.4), and 
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Figure 11: (Left) F — gJ 2 F$ + ^logfc is plotted against l/N 2 . The results are nicely fitted to 
straight lines, which enables us to extract the genus one contribution reliably. (Right) The genus 
one contribution extracted from the left panel is plotted against A. The solid line represents the 
genus one contribution to -Ffhm > whereas the dotted line represents -FpHM — c i> wnere c i is given 
by (5.15). The dashed line represents the weak coupling behavior given by (5.13) with the -\\ogk 
term being subtracted. 



subtract also the term — glog/c that appears in the weak coupling result (5.13). Then we 
plot the result after these subtractions against l/N 2 in fig. 11 (Left), which can be nicely 
fitted to straight lines. The intercepts give the values of the genus one contribution for 
each A, which are plotted against A in fig. 11 (Right). We find that the result is in good 
agreement with the genus one contribution of -Ffhm after making a constant shift by —c\ 
determined as (5.15). 

5.2 Finite N effects 

Let us see how well the FHM result with the corrections (5.18) does at finite N. In fig. 12 
we plot F — FfhMj be., the discrepancies between our result and the FHM result for N = 2 
and N = 10 against k. At k > 1, the N = 2 data and the N = 10 data are on top of each 
other as anticipated from fig. 8. In this regime, the discrepancies are in good agreement 
with the corrections (5.18) identified in section 5.1. 

It is interesting to see what happens if we go to smaller k region in fig. 12 although non- 
integer k is not physical in the original gauge theory. Firstly we start to see some difference 
between N = 2 and N = 10, which implies that there is some N dependence which is not 
captured by the FHM result in this regime. We speculate that the N dependence is due to 
the membrane instanton effect [31, 57], which behaves as 0(e~ 7Ty " 2kN ), since the worldsheet 
instanton effect is negligible in this regime. Secondly, we find that the discrepancy between 
our result and the FHM result no longer agrees with (5.18) including corrections up to 
genus two. This is understandable since the higher genus terms become non- negligible as 
one goes to smaller k (larger g s ). 

On the other hand, the free energy at small k is calculated by the Fermi gas approach 
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Figure 12: (Left) The discrepancy of the free energy between our data and the FHM result is 
plotted against k, and compared with the analytical results around the planar limit for N = 2 
and N = 10 . The dashed and the solid lines represent the correction (5.18) up to genus one 
term and up to genus two term, respectively. The dotted line represents the sum of the constant 
map contributions at all genus (5.19). (Right) The same as the left panel except that our data 
is compared with the result obtained by the Fermi gas approach. The dotted and the solid lines 
represent the analytical results (2.12) with A(k) given by (2.11) and (2.13), respectively. The dotted 
line represents again the sum of the constant map contributions at all genus (5.19). 

as (2.12). We find that our data for N = 10 interpolate nicely the behavior (2.12) at 
small k and the behavior (5.18) at large k. This also supports our speculation that the 
difference between the N = 2 and N = 10 data at small k is due to the membrane instanton 
effect, which is neglected in the Fermi gas approach. Note, in particular, that the Fermi 
gas approach yields correction to the FHM result in odd powers of k, whereas the genus 
expansion yields even powers of \jk. Our data nicely interpolate the two asymptotic 
behaviors, which are smoothly connected around k ~ 1. 

Finally let us consider the sum of the constant map contributions at all genus (5.17), 
which is given by 13 

Fconst = "|^ 2 " \ fog k + 1 log | + 2C'(-1) 

1 f°° 1/31 3 \ , 

-o ^i^i l r-rr- • (5-19) 

3 Jo e ~~ 1 \ x x xsmh x J 

In fig. 12 we find that this function agrees well with the discrepancy between our result 
and the FHM result for the whole range of k investigated, including k < 1. Since the Fermi 
gas result (2.12) also gives accurate description there, it is natural to guess that they are 
actually the same. Indeed, as we show in appendix C, the sum of the constant map (5.19) 
can be expanded around k = as 

_ 2C(3) 1 \ ^ ("I)" p R 2n-2.2n-l 

r const — 51 7T lo S 2 + / , ,r> \ i -P2n.-P2n-27T K 

n=i y ' 



An assumption is needed to obtain (5.19) and (5.20). For details, see appendix C. 
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1 n 2C(3) k TT 2 k 3 7T 4 fc 5 , rort , 

which exactly reproduces the result of the Fermi gas approach (2.12) to the A; 3 term. 14 
Remarkably, the constant map contribution (5.19) and the expansion (5.20) are equivalent 
at any k. Therefore, we expect that the result (5.20) is the all-order form of A(k) — \ log 2 
in the Fermi gas approach. In other words, we expect that the expansions of the free energy 
around k = oo (the constant map contribution) and k = (the Fermi gas approach) give 
the same answer after taking sums to all orders. In this sense the free energy around the 
planar and M-theory limits are smoothly connected with each other. 

6. Summary and discussions 

In this paper we have established a simple numerical method for studying the ABJM theory 
on a three sphere for arbitrary rank N at arbitrary Chern-Simons level k. The crucial point 
is that we are able to rewrite the ABJM matrix model, which is obtained after applying the 
localization technique, in such a way that the integrand becomes positive definite. By using 
this method, we have confirmed from first principles that the free energy in the M-theory 
limit grows proportionally to TV 3 / 2 as predicted from the eleven-dimensional supergravity. 
We have also found that the FHM formula with the additional terms describes the free 
energy of the ABJM theory in the type IIA superstring and M-theory regimes. Analytic 
form of the additional terms is discussed in detail, and beautiful agreement between planar 
and M-theory regions is found. These additional terms become important when we consider 
the quantum string effect in the AdS/CFT duality. 

There are many issues worth being addressed by using our method. Most importantly 
from the conceptual point of view, we can use the free energy obtained for finite TV and 
finite k to test the AdS/CFT duality at the quantum string level. At the tree level, or 
equivalently in the planar limit, there is strong evidence that the gauge theory correctly 
describes the string ol effect. For example, the DO-brane quantum mechanics reproduces 
the a 1 correction to the black 0-brane solution in type IIA superstring theory [60]. However, 
no definite conclusion is obtained for quantum string corrections so far. In fact, as pointed 
out in ref. [33], the FHM formula does not seem to agree with a prediction from the 
string theory side [49]. This disagreement is not solved even if we take into account of the 
corrections found in this paper. Some of the possible solutions to this puzzle includes: (i) 
one has to consider some different gauge theory such as SU(N)k x SU(N)-k theory, which 
gives different l/N corrections, (ii) one has to refine the argument on the string theory 
side, and (iii) the AdS/CFT does not hold at the quantum string level. In particular, the 
scenario (i) can be tested straightforwardly by extending our method. 

We consider it equally important to study quantum M-theory. While there is very 
little knowledge on it so far, we may hope to get some insight through intensive numerical 

14 Note that the Fermi gas calculation has been done only to the fc 3 term, and the higher order terms are 
just an educated guess. Calculations at higher orders would be an interesting future direction. 
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studies of the ABJM theory. In fact similar attempts have been made recently [61, 62] 
using the BFSS matrix theory [10]. Numerical studies suggest that the prediction from 
the dual string theory for the scaling dimension of a certain class of operators continues to 
hold in the M-theory region. Similar or possibly more striking features of M-theory may 
show up by studying the ABJM theory numerically. 

While we have focused on the free energy as the most fundamental quantity in the 
ABJM theory, our method can be used to calculate the expectation values of BPS operators. 
For instance, it is possible to calculate the expectation value of the circular Wilson loop for 
various representations. They are conjectured to be related to the string worldsheet area 
and the D-brane world-volume in the type IIA superstring region, respectively. It would 
be interesting to test this conjecture and to see the stringy corrections. 

Our method can be also applied to other theories, which have recently attracted much 
attention in connection to the F-theorem and the entanglement entropy. For example, 
one can study the necklace- type quiver discussed in ref. [63]. We can also study other 
gauge groups such as O(N) and USp(N) studied in ref. [64, 65]. Detailed studies of these 
theories outside the planar limit, in particular, would be very interesting. For instance, the 
so-called orbifold equivalence, which is usually proven to hold only in the planar limit, can 
hold outside the planar limit in these models [66, 67]. Note also that the ABJM matrix 
model is related to the lens space matrix model, which appears in the context of the 
topological string theory. It is therefore conceivable that there might be some applications 
to the topological string theory as well. 

We hope that the results of this work are convincing enough to show the power of the 
numerical approach, and that many more applications other than those listed above would 
reveal themselves as we proceed further. 
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A. Basics and details of the Monte Carlo simulation 



In this section we explain the basics and details of the Monte Carlo simulation 15 . Let us 
consider the action 

/n,v, tanh 2 ((xi -Xj)/2k)\ 

SW ^' > = -'°H n.2cI fe /2) - )• < A1 > 

(Below we abbreviate S(N, k; x) as S(x).) Let 0(x) be an "observable", which is a function 
of {xi}. In general, it is difficult to calculate the expectation value of O defined by 

fd N xO(x)e' s ^ IK . 

< 0)= s^J-m ■ (A - 2) 

A brute force integration is not practical unless the number of variables N is very small 
such as N < 5. Monte Carlo simulation is a practical tool, which enables this calculation 
even for large N. 

In Monte Carlo simulation, a series of configurations {xj} 

K (0) } {xf ] } - (A.3) 

is generated in such a way that the probability with which a configuration {x{\ appears 
approaches e~ s ^/Z as the number of configurations increases. More precisely, we re- 
quire that the probability Wk{{xi}) with which a configuration {xi} appears at M-th step 
converge to e~ s ^/Z as 

£ -S(x) 

lim w M ({xi}) = —=— . (A.4) 

M-»oo Z 

Then the expectation value can be obtained by simply taking an average over the config- 
urations {a^- n ' ) } as 

1 M 

(0(x)}= lim -^O^). (A.5) 

n=l 

This can be achieved by generating the series with a transition probability P({x^} — > 
{a^ n+1 ^}), which (neglecting a few technical details) satisfies following conditions. 

• Markov chain. — The transition probability from {x^} to {xj n+1 ^} does not depend 
on previous configurations {x[^} (I < n). 

• Ergodicity. — For any pair of configurations {x} and {x 1 }, there is nonzero transition 
probability within finite steps. 

• Aperiodicity. — The probability from {x^ to {x^ is always nonzero. 



15 There are many good references on Monte Carlo methods. See, e.g., ref. 
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• Positive state. — All configurations have finite mean recurrence time 16 . 

• Detailed balance. — The following equality should hold for arbitrary pairs of config- 
urations {x} and {x'}. 



e - s( -^P{x -> x 1 ) = e- s ^P{x' -> x) . (A.6) 



There are many ways to satisfy these conditions. In this work we use the Hybrid 
Monte Carlo (HMC) method [59]. We introduce fictitious momentum variables pi (i = 
1, 2, ■ ■ ■ , N), which are conjugate to Xj, and consider a Hamiltonian 

h = J2^ + s{x) - (A - 7) 

i 

Starting with an initial configuration {xf^}, we generate a series of configurations {x^} (n - 
1, 2, ■ ■ ■ ) by repeating the following steps: 

• Generate vT~^ randomly, with a probability weight e ( Pl ) ^ 2 . 



• Starting with a configuration {xi,p{\ = {x\ n l \p^ get a new configuration 
{x^p'j} by the "molecular dynamics" explained below. 

• "Accept" the new configuration {x^p'j} (i.e. take {x^} = {x^}) with a probability 
min{l,e^~^ }, where H and H' are the value of the Hamiltonian calculated with 
{xi,pi} and {x'iip'j}, respectively. When the new configuration is rejected, we keep 
an old configuration, so that {x^} = {xi} = {xj™" 1 -*}. 

The "molecular dynamics" is defined as follows. First we introduce a fictitious time r 
and consider the time evolution according to the Hamilton equations 

dxi dpi dH dpi dH dS 

dr dr dpi ^ % ' dr dxi dxi 

If we solve the Hamilton equations exactly, the Hamiltonian is conserved. In practice, we 
solve them approximately by discretizing the differential equations, so the Hamiltonian is 
not conserved exactly. We denote the time step as At and the number of steps as N T . 
Then x\ = xi(N t At) and p\ = pi(N T Ar) are calculated by using the following "leap-frog 
method", starting with Xj(0) = x, and Pi(0) = p%. 



. Xi(^) = Xi (0) + Pi (0) • % 



• Pi (Ar) =p<(0) - 



dS 



■At 



2 



• x % (|Ar) = Xi (^) + Pi (At) • At 

16 If P/™' is the probability to get from {x} to {x} in n-steps of the Markov chain, without reaching this 
configuration at any intermediate step, then the mean recurrence time of {x} is defined by r, = YL^=i n ^ii™' ) • 
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Pi (2Ar)= Pi (AT)- w 



s -Ar 

At 



Xi ((N T - 1/2) At) = Xi ((N T - 3/2) At) + Pi ((N T - l)Ar) • Ar 



p i {N T AT)=p l {{N T -l)AT)- £ 



• Ar 

=(JV t -1/2)At 

Ar 



Xi {N T A T ) = Xi ((N T - 1/2) At) + Pi((N t )At) ■ %f 



Note that the leap-frog method is designed so that the reversibility is satisfied. Namely, 
by starting with the final configuration {x^} and {p[} and reversing the time, the initial 
configuration {xj} and {pi} is reproduced. As a result, the detailed balance condition is 
satisfied. 17 

In the simulation, N T and Ar should be chosen so that a good approximation is 
achieved with fewer configurations. For that purpose, (i) the change at each transition 
should be sufficiently large, and (ii) the acceptance rate should be large. The first condition 
is achieved by taking N t At sufficiently large. However, if we fix Ar and take larger iV T , 
the Hamiltonian is not conserved at all, and the new configurations are hardly accepted. 
Therefore one has to take Ar smaller so that the conservation of the Hamiltonian becomes 
better. In actual simulations (at N = 20 and k = 5, for example), we took Ar ~ 0.1 and 
N T ~ 200, so that the acceptance rate is around 0.8. As an initial configuration, we choose 
x\°^ to be a random number in [—0.5,0.5]. 

In Monte Carlo simulation, configurations with larger path-integral weight ("important 
samples") appear more often. For this reason it is called also the importance sampling. 
Since the region of configuration space, which gives dominant contribution to the path 
integral is typically quite limited, a good approximation can be achieved with a rather 
small number of important samples. This should be contrasted to the usual brute force 
integration, in which most of the CPU time is wasted for calculating the integrand for 
unimportant configurations. 



17 For simplicity, let us assume H > H' . (The argument for the case with H < H' is the same.) Because 
of the reversibility, the transition probabilities are 

P{{ Xi ] -> {x[}) = e- p2/2 /^ x min{l, e H - H '} (A.9) 



= e 



-p 2 /2 



/Vn (A.10) 



and 



P({x'i} -> {Xi}) = e- p ' /2 /v^ x min{l, e 11 '-"} (A.ll) 

= e-^'l-W/^. (A.12) 

Therefore, 

e- s(x) P({x t } -> {x[}) = e- s{x) x e- p2/2 /V^ (A.13) 

= e - s (*') P ({^} -»■ { Xi }) . (A.14) 



- 25 - 



In Monte Carlo simulation, as we have described above, configurations are generated 
with the probability e~ /Z. Therefore, the Monte Carlo method works only if the path- 
integral weight e~ s is real and positive. If the measure e~ s is not real and positive, the 
model is said to suffer from the sign problem or the phase problem; here "sign" and "phase" 
mean the negative sign and the complex phase of the integration weight. In the original 
form of the ABJM matrix model (1.2), the partition function is given by an integration 
of a complex function. Therefore it suffers from the sign problem, and the Monte Carlo 
method is not applicable. 



B. Derivation of the sign-problem-free form of the ABJM matrix model 

In this section we derive the sign-problem-free form of the ABJM matrix model, which was 
used in ref. [51, 34, 35] for a different purpose. Let us start with the ABJM matrix model 
(1.2). We are going to use the Cauchy identity 18 

Ui<j( U i -Uj)(Vi -Vj) 



Here a runs through all permutations. By setting U{ = e Mi ,t>j = e"*, it becomes 
n^-e^-e^') = £ ( _ ir jj 1 



(B.l) 



e m _|_ e Mo 



From this, we obtain 



rii<j 


2sinh(^ l / J )] 


[2sinh(^)] 




"2cosh(^)" 





= E<-i>TI 



2 cosh 



Therefore, the partition function can be written as 



Z(N,k) 

J_ Vf-ir +CT ' 



d N ix d N u 



n 



i 



W V [2 cosh (Z=F£L) 

_ J_ vr-ir f dNfX dNv TT 

~ ) J (2^(2^11 [ 2cosh( «) 

By using the formula 



2 cosh 



2 cosh 



exp 



— = -[ 

2 cosh p it J 



alx 



e T 
2 cosh x ' 



we obtain 



E(-D'n 



r 2c osh (i^)l f2cosh (e^sa) 



3 See the appendix of ref. [51] for the proof of this identity. 
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= EM)^ / dVy 

(J 



exp 



exp 



Y\i 2 cosh • 2 cosh 



Hi 2 cosh Xj • 2 cosh j/j 



Therefore, the partition function becomes 19 



2 cosh Xj • 2 cosh 



d^/i d^ 



i i ik N 

exp [~ ~ "O 2 * + ~ E^ iyi ~ + ^ E(^ 2 ~ K < 

i i i=l 



(2tt) w (2tt) w 

L vc-iv 7 — / d^xd^w — 

TV! ' ir 2N J 1 1 ■> .•)„>,!,„ / ro^i.v 



exp 



d"ji__cFv 
2 cosh Xi • 2 cosh y 4 7 (2^)^ (2vr) 



A? 



i=l 



A? 



i=l 



. AT 

exp [- E + yi ) 2 ~ + y<ja)) 2 ) 



i=l 



{ ik" 2 ik" 2 2*^ , 
ex P Li^ E ^ " 4^ E ^ " ^ E " 

i=l i=l i=l 



d^V 

, 2 cosh Xi ■ 2 cosh Vi J (2ir) N (2ir) N 

N 



2 cosh yi 



1 Vf-ir - ! d N xd N v — 

<J 1 



■|^E£l^(S/i-2/a(i)) 



a 2 cosh (f). 2 cosh (f-Jtfa) 



(B.4) 



19 In the Fermi gas approach [35], the integrand is identified with a partition function for the ideal Fermi 
gas given by 



1 r N 

z(N,k) = —^(-wj ^ Yi P ( Xi ,x a{i) ), 



(B.2) 



where p(x\,X2) is interpreted as the one-particle density matrix 

11 1 

P(Xl, x 2 ) 



2vk (2cosh^) 1/2 (2coshf ) 1/2 2cosh(^; 



(B.3) 
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We use the Cauchy identity again: 

f2 S inh(^ 



£(-irII 



, 2 rush (^-) Ui,j 
Thus we arrive at the final expression 



2 cosh 



2fc ) 



^ri tanh2 ( : 

i<j v 



2A; 



Z(Nk)- 1 /" ^ n^tanh^^ 
Z(7V ' k) -WW.J (2^F n,2cosh(f) ' (B - 5) 

which does not have a sign problem. 

C. The relation between the constant map and the Fermi gas result 

In this appendix we show the correspondence between the constant map contribution and 
the Fermi gas result A(k) — \ log 2, which is derived by the large-A; and small-fc expansions, 
respectively. 

As we mentioned earlier, the constant map contribution F const is given by 

oo 

Fconst=E^ 2 ^ott, (CI) 
9=0 

where the coefficients F^ st are 

p(0) _ C(3) = of(_-\\ + llntrJL F (3> 2 ) _ 19 B 2 gB 2 g-2 (r g s 

const 2 . const ^ 6 S 2fc ' const (4 5 ) (2g - 2) (2 5 - 2) ! k ; 

In order to evaluate the summation more easily, we use the integral representation of the 
Bernoulli number, 

»oo 2g-l 



B 2 g 



roo -Zg-L 

{ - l)9 ~ H9 j Q ^—i dx (9 = 1,2,---)- (C3) 
By using this representation, we obtain 

oo oo j-> poo 2o— 3 

Z^9s ^const L\ 9 S 4 (2 5 )(2 5 -2)! 7o e 2 ™-l 

= „-2 / Y(-1) B , 9 ^ 9sXf 9 • (C.4) 

3 Jo e 2 «-l^ 1 {2g){2g-2)\ K ys 1 K ' 

g—z 

This summation is easily performed by using a formula 

Note that the series converges only for \z\ = \g s x\ < it. However, since the result of the 
summation in the right hand side is well-defined even for \g s x\ > n, we analytically continue 
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g s x to the whole region including \g s x\ > tt and assume that this does not affect the result 
of the integration 20 . 

By substituting z = g s x in the formula, the summation is rewritten as 

V s - 2o ,„M * 2 f°° , I 4T 2 , 12ir 2 i 2 \ ,„ , 

By changing the variable as t = we obtain a simpler form, 

= -\ f * ^ (a - f - ^) • (c.7, 

g—2 

Although each term of the integrand is divergent at t ~ 0, this is canceled with each other, 
and therefore the integral gives a finite value. In order to make our analysis easier, we will 
apply the zeta-function regularization to the integral. 
For later convenience, we decompose the integral as 

f>?-^±, = f * ^ (-£ + i) + if * • (cs) 

Note that the first factor in the second term is the generating function of the Bernoulli 
number 



e 

n=0 



Although this series also converges only for \x\ < 2ir, we analytically continue it to the 
whole region and assume that this does not affect the result. Then by using the formula 



e x _ I jL^ sm tf x jL^ 

m=l m=l 

the integral rewritten as 

00 poo / -I 1 \ /I 00 R 00 foo 

s =2 m=l u v 7 ra=0 m=l ,/u 

(C.ll) 

The first integral is easily performed by using 

f CO /'CO 

/ tft e" 5 **" 1 = -7 -logs, / dt e" st t" 3 = --s 2 (-3 + 2 7 + 21ogs) , 
where 7 is the Euler-Mascheroni constant. We obtain 



00 „ c 



Even if this assumption is not valid, the discrepancy at small k is expected to be c at most 
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= — £ n 2 (-3 + 2 7 + 2 log (km)) + - £ (-7 - log (mk)) 

m=l m=l 

= ^ [(-3 + 2 7 + 2 log fc)C(-2) - 2C'(-2)] + 1 [(-7 - log fc)C(O) + C'(0)] 
= -y C'("2) + I [(7 + log fc) - log (2tt)] . (C.12) 
Next we evaluate the second integral 

OO „oo 



^ PC 

> ml 

^ Jo 



dt t~^ n e 



2+n-2mt 



For n = 

By using the formula 

poo 



poo 

/ dt e^tT 2 = s(-l + 7 + logs) , (C.13) 
Jo 



we obtain 



00 „ c 
> m I 



dt t- 2 e~ 2mt = -2C'(-2) . (C.14) 



• For n = 1 



00 />oo 1 

E ffl dtt- 1 e- 2mt = -(7 + log2) + C , (-l). (C.15) 

m=l ^° 



• For n > 2 

By using the formula 



/•OO 1 

/ e" 5 *^ 1 = T(A)— (Re(X) > 0) , (C.16) 
jo s 



the integral becomes 



'OO 



£>/ d^-^ e - 2 -* = £^-ilc(n-2). (C.17) 

m=l 7 ° ^ 
Thus the constant map contribution is rewritten as 

/■/Q") 1 1 7T fc 2 1 

Fconst = -f^k 2 - - Jog fc + - log - + 2C'(-1) - yC'(-2) + - (7 + log k - log (2vr)) 



4 
+ k 



-2B oC '(-2) + Bxife (1( 7 + log 2) + C '(-l)) + £ ^ fe2n£ |tA(2n - 2) 
( 8^ + ^T^) + ( 2(1 + 2B ^'^ + ^ + 2Bl ^ + ^(-! + 5 i) lQ g 2 ) 
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o R 
--B^-2) + g ^-^_ C( 2» - 2)*— . (C.18) 

Since B„ = 1, B t = - >, f(-2) = -@ and f(2n) = Hr'^fn!., we obtain 

11=1 ^ ' 

1, „ 2C(3) k A 3 ir 4 t 5 
= -2 log2 + l^-l2-4350 + 907200 + ---' (C - 20) 

which is same as A{k) — \ log 2 derived by the Fermi gas approach [35] up to the order of 
0(/c 5 ). Therefore, we expect that this is the all-order form in the Fermi gas picture if we 
calculate higher order of k. 

Thus we conclude that the constant map contribution and the term A(k) — \ log 2 in 
the Fermi gas result are the series expansions of the integral representation (C.7) around 
k = oo and k = 0, respectively, with the radius of convergence being finite. In other words, 
the two expansions are smoothly connected with each other by analytic continuation. 
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